function DJ = s_grad_2(model,fid,kbasis,cgp)

% S_GRAD Compute gradient of error function
%
% MODEL is the parameter array
% FID is the FID array
% KBASIS is a basis array defined by make_kbasis().
% CGP is the reconstruction parameter structure.

if isfield(fid,'N') 
    NF = fid.N;
else
    NF = length(fid.sig);
end
if isfield(model,'N') 
    NM = model.N; 
else
    NM = length(model.pvec)/3; 
end

sigdiff1 = fid.sig1(1:NF) - fid.est1(1:NF);
sigdiff2 = fid.sig2(1:NF) - fid.est2(1:NF);
%figure(16)
%plot(real(sigdiff))
%drawnow
wfr = model.wfr;
wfi = model.wfi;
ix = model.ix;
iy = model.iy;
alphaM = cgp.alphaM;
alphaf = cgp.alphaf;
if (isreal(alphaf))
    alphaf = (1+1i)*alphaf;
end
p_image = zeros(model.reso);
p_image1 = zeros(model.reso);p_image2 = zeros(model.reso);
% p_grad = zeros(model.reso);
p_image1(ix+model.reso*(iy-1)) = model.pvec(1:NM);
p_image2(ix+model.reso*(iy-1)) = model.pvec(NM+1:2*NM);
% M_grad1 = zeros(NM,1);
% M_grad2 = zeros(NM,1);
f_grad1 = zeros(NM,1,'single');
f_grad2 = zeros(NM,1,'single');
CC_conj1 = conj(-sigdiff1);
CC_conj2 = conj(-sigdiff2);
L = [ 0 1 0; 1 -4 1; 0 1 0];
M1 = transpose(model.pvec(1:NM));
M2 = transpose(model.pvec(NM+1:2*NM));
wf = weight_cexp(model.pvec(2*NM+1:3*NM).',wfr+1i*wfi);
% wf = model.pvec(2*NM+1:3*NM).';
wf_n = ones(NM,1,'single');

if isstruct(kbasis)
    % use separable calculation
    M_grad1 = (kbasis.ekxx(1,ix).*kbasis.ekyy(1,iy)).'*CC_conj1(1);
    M_grad2 = (kbasis.ekxx(1,ix).*kbasis.ekyy(1,iy)).'*CC_conj2(1);
else
    M_grad1 = kbasis(:,1)*CC_conj1(1);
    M_grad2 = kbasis(:,1)*CC_conj2(1);
end

if (alphaM ~= 0)
    p_grad1 = conv2(conv2(conj(p_image1),L,'valid'),L);
    M_grad1 = M_grad1 + alphaM*p_grad1(ix+model.reso*(iy-1));
    p_grad2 = conv2(conv2(conj(p_image2),L,'valid'),L);
    M_grad2 = M_grad2 + alphaM*p_grad2(ix+model.reso*(iy-1));
    if (1==0)
        figure(16)
        subplot(1,2,1)
        imagesc(real(p_image1))
        subplot(1,2,2)
        imagesc(real(p_grad1))
        drawnow
        pause;
    end
end


wf_n = wf_n.*wf;
for ii = 2:NF   %#ok<FORPF>
    
     if isstruct(kbasis)
         % use separable calculation
         Cekxx1 = kbasis.ekxx(ii,:)*CC_conj1(ii);
         temp1 = wf_n.*(Cekxx1(ix).*kbasis.ekyy(ii,iy)).';
         Cekxx2 = kbasis.ekxx(ii,:)*CC_conj2(ii);
         temp2 = wf_n.*(Cekxx2(ix).*kbasis.ekyy(ii,iy)).';
         M_grad1 = M_grad1 + temp1;
         M_grad2 = M_grad2 + temp2;
         f_grad1 = f_grad1 +(ii-1).*temp1;
         f_grad2 = f_grad2 +(ii-1).*temp2;
     else
         temp0 = wf_n.*kbasis(:,ii);
         temp1 = temp0*CC_conj1(ii);
         temp2 = temp0*CC_conj2(ii);
         M_grad1 = M_grad1 + temp1;
         M_grad2 = M_grad2 + temp2;
         f_grad1 = f_grad1 +(ii-1).*temp1;
         f_grad2 = f_grad2 +(ii-1).*temp2;
     end

end

f_grad = M1.*f_grad1 + M2.*f_grad2;

if (alphaf ~= 0)
    p_image(ix+model.reso*(iy-1)) = model.pvec(2*NM+1:3*NM);
    p_grad = conv2(conv2(conj(p_image),L,'valid'),L);
    p_grad = wfr*real(p_grad) + wfi*imag(p_grad);
    pg = p_grad(ix+model.reso*(iy-1));
    f_grad = f_grad + wfr*real(alphaf)*real(pg)+1i*wfi*imag(alphaf)*imag(pg);
end

DJ = [2*conj(M_grad1);2*conj(M_grad2);2*wfr*real(f_grad)-1i*2*wfi*imag(f_grad)].';
